%% My First Flow Solver: Gravity Column
% In this example, we introduce a simple pressure solver and use it to
% solve the single-phase pressure equation
%
% $$\nabla\cdot v = q, \qquad
% v=\textbf{--}\frac{K}{\mu} \bigl[\nabla p+\rho g\nabla z\bigr],$$
%
% within the domain [0,1]x[0,1]x[0,30] using a Cartesian grid with
% homogeneous isotropic permeability of 100 mD. The fluid has density 1000
% kg/m^3 and viscosity 1 cP and the pressure is 100 bar at the top of the
% structure.
%
% The purpose of the example is to show the basic steps for setting up,
% solving, and visualizing a flow problem. More details on the grid
% structure, the structure used to hold the solutions, and so on, are given
% in the <simpleBC.html basic flow-solver tutorial>.
try
require incomp
catch %#ok<CTCH>
mrstModule add incomp
end
%% Define the model
% To set up a model, we need: a grid, rock properties (permeability), a
% fluid object with density and viscosity, and boundary conditions.
gravity reset on
G = cartGrid([1, 1, 30], [1, 1, 30]);
G = computeGeometry(G);
rock.perm = repmat(0.1*darcy(), [G.cells.num, 1]);
fluid = initSingleFluid('mu' , 1*centi*poise, ...
'rho', 1014*kilogram/meter^3);
bc = pside([], G, 'TOP', 100.*barsa());
%% Assemble and solve the linear system
% To solve the flow problem, we use the standard two-point
% flux-approximation method (TPFA), which for a Cartesian grid is the same
% as a classical seven-point finite-difference scheme for Poisson's
% equation. This is done in two steps: first we compute the
% transmissibilities and then we assemble and solve the corresponding
% discrete system.
T = computeTrans(G, rock);
sol = incompTPFA(initResSol(G, 0.0), G, T, fluid, 'bc', bc);
%% Plot the face pressures
clf
plotFaces(G, 1:G.faces.num, convertTo(sol.facePressure, barsa()));
set(gca, 'ZDir', 'reverse'), title('Pressure [bar]')
view(3), colorbar
set(gca,'DataAspect',[1 1 10]);
%%
displayEndOfDemoMessage(mfilename)The rest of the example will walk through this example line for line, displaying the PRST examples.
In [1]:
%load_ext line_profiler
In [2]:
# For Python 3 compatibility
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import numpy as np
# Enable MayaVi to work non-blocking in notebook
#import matplotlib
#matplotlib.use('Qt4Agg')
#matplotlib.interactive(True)
import prst
import prst.incomp as incomp
import prst.gridprocessing as gridprocessing
import prst.utils as utils
import prst.params as params
import prst.solvers as solvers
import prst.plotting as plotting
from prst.utils.units import *
In [3]:
help(prst)
%% Define the model
% To set up a model, we need: a grid, rock properties (permeability), a
% fluid object with density and viscosity, and boundary conditions.
gravity reset on
In [4]:
# PRST uses a single module-level variable to set the gravity
prst.gravity_reset()
print("Gravity vector:", prst.gravity)
G = cartGrid([1, 1, 30], [1, 1, 30]);
G = computeGeometry(G);
rock.perm = repmat(0.1*darcy(), [G.cells.num, 1]);
fluid = initSingleFluid('mu' , 1*centi*poise, ...
'rho', 1014*kilogram/meter^3);
In [5]:
G = gridprocessing.cartGrid([1, 1, 30], [1, 1, 30])
print(G)
Computing the geometry modifies the original grid in place, no assignment is necessary. Note that specifying the module is necessary in PRST. This avoids polluting the global namespace, at the cost of more verbosity.
In [6]:
gridprocessing.computeGeometry(G)
Out[6]:
In [7]:
rock = params.rock.Rock(G, perm=0.1*darcy, poro=1)
In [8]:
fluid = incomp.fluid.SingleFluid(viscosity=1*centi*poise, density=1014*kilogram/meter**3)
print(fluid)
Set up boundary conditions. Constant pressure at top face.
In [9]:
bc = params.wells_and_bc.BoundaryCondition()
bc.addPressureSide(G, "top", 100*bar)
Out[9]:
In [10]:
print(bc)
%% Assemble and solve the linear system
% To solve the flow problem, we use the standard two-point
% flux-approximation method (TPFA), which for a Cartesian grid is the same
% as a classical seven-point finite-difference scheme for Poisson's
% equation. This is done in two steps: first we compute the
% transmissibilities and then we assemble and solve the corresponding
% discrete system.
T = computeTrans(G, rock);
sol = incompTPFA(initResSol(G, 0.0), G, T, fluid, 'bc', bc);
In [11]:
T = solvers.computeTrans(G, rock)
resSol = solvers.initResSol(G, p0=0.0)
sol = incomp.incompTPFA(resSol, G, T, fluid, bc=bc)
In [12]:
sol
Out[12]:
%% Plot the face pressures
clf
plotFaces(G, 1:G.faces.num, convertTo(sol.facePressure, barsa()));
set(gca, 'ZDir', 'reverse'), title('Pressure [bar]')
view(3), colorbar
set(gca,'DataAspect',[1 1 10]);
In [13]:
import prst.plotting as plotting
from prst.utils.units import convert
from mayavi import mlab
plotting.plotGrid(G, cell_data=convert(sol.pressure, to=bar),
colorbar=True, mlab_show=False)
mlab.axes()
mlab.show()
Out[13]: